;;----------------------------------------------------------------------
;; plot zonal distribution of ion production 
;;----------------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$KAILIB_ROOT/lib_radon.ncl"

begin

  plotname = "fig_emis" 
  plotform = "eps" 
  plotfile = plotname+"."+plotform 
  plotpdf  = plotname+".pdf"

  StringFontHeightF  = 0.030 
  cnLineLabelFontHeightF = 0.025 
  tmXBLabelFontHeightF = 0.022
  tmYLLabelFontHeightF = 0.022
  tiMainFontHeightF  = 0.035 
  tiXAxisFontHeightF = 0.035
  tiYAxisFontHeightF = 0.035
  lbLabelFontHeightF = 0.030 

  gsnStringFont = "helvetica-bold"
  gsnStringFont = "helvetica"

  np = 6

  str = (/"","","","","","","","",""/) 

  tiYAxisString = str
  tiYAxisString(0) = "Emission rate (atom m~S~2~N~ s~S~-1~N~) " 
  tiYAxisString(3) = "Emission rate (atom m~S~2~N~ s~S~-1~N~) " 
  tiYAxisString(6) = "Emission rate (atom m~S~2~N~ s~S~-1~N~) " 

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/) 

  abcd = abcd + ") " 

  gsnLeftString = str

  stra = "WCRP1995 " 
  strb = "SW1998 scaled " 
  strc = "Merged " 
 
  strd = "Data source of the merged map" 
  stre = "MERGE land " 
  strf = "MERGE ocean " 
 
  gsnLeftString(0) = abcd(0) + stra
  gsnLeftString(1) = abcd(1) + strb
  gsnLeftString(2) = abcd(2) + strc

  gsnLeftString(3) = abcd(3) + strd 
  gsnLeftString(4) = abcd(4) + stre
  gsnLeftString(5) = abcd(5) + strf

  gsnLeftString(6) = abcd(6) + stra 
  gsnLeftString(7) = abcd(7) + strb
  gsnLeftString(8) = abcd(8) + strc 
 
  gsnCenterString = str

  strx = "" ;;"         IPR ( cm~S~-3~N~s~S~-1~N~ ) "
  stra = strx ;;+ "WCRP1995 " 
  strb = strx ;;+ "SW1998 scaled " 
  strc = strx ;;+ "Merged " 

  gsnCenterString(0) = stra
  gsnCenterString(1) = strb
  gsnCenterString(2) = strc 

  gsnCenterString(3) = stra 
  gsnCenterString(4) = strb
  gsnCenterString(5) = strc 

  gsnCenterString(6) = stra 
  gsnCenterString(7) = strb
  gsnCenterString(8) = strc 

  gsnRightString = str
 
  stra = "ANN" 
  strb = "ANN" 
  strc = "ANN"

  gsnRightString(0) = stra
  gsnRightString(1) = stra
  gsnRightString(2) = stra 

  gsnRightString(3) = strb 
  gsnRightString(4) = strb
  gsnRightString(5) = strb 

  gsnRightString(6) = strc 
  gsnRightString(7) = strc
  gsnRightString(8) = strc 

;;----------------------------------------------------------------------
;; load files 
;;----------------------------------------------------------------------

  overland = True
 
  path = "$KAILIB_ROOT/ionprod/"
  pata = "/pf/m/m222044/rnpaper/"


  run1 = "rnemis_E0" 
  run2 = "rnemis_E2" 
  run3 = "rnemis_E3" 
  runa = "rnemis.E91.T63.clim.nc" 
  runb = "rnemis.E92.T63.clim.nc" 
  runc = "rnemis.E93.T63.clim.nc" 
  rund = "rnemis.E94.T63.clim.nc" 
  rune = "rnemis.E95.T63.clim.nc" 

  filenama = run1+"_ANN.nc"   
  filenamf = run2+"_ANN.nc"   
  filenami = run3+"_ANN.nc"  

  filenamj = runa+".nc"  
  filenamk = runb+".nc"  
  filenaml = runc+".nc"  
  filenamm = rund+".nc"  
  filenamn = rune+".nc"  

  filenamd = path+"mbase_T63L31.nc" 

  fila = addfile(filenama,"r") 
  filf = addfile(filenamf,"r") 
  fili = addfile(filenami,"r") 

  fild = addfile(filenamd,"r")

  filj = addfile(filenamj,"r")
  filk = addfile(filenamk,"r")
  fill = addfile(filenaml,"r")
  film = addfile(filenamm,"r")
  filn = addfile(filenamn,"r")

;;----------------------------------------------------------------------
;; read data 
;;----------------------------------------------------------------------

  rna = fila->rnemis
  rnf = filf->rnemis
  rni = fili->rnemis

  rnj = filj->rnemis
  rnk = filk->rnemis
  rnl = fill->rnemis
  rnm = film->rnemis
  rnn = filn->rnemis

  eps = 1.E-05

  rnj = where(ismissing(rnj),-1,1.-eps) 
  rnk = where(ismissing(rnk),-1,2.-eps) 
  rnl = where(ismissing(rnl),-1,3.-eps) 
  rnm = where(ismissing(rnm),-1,4.-eps) 
  rnn = where(ismissing(rnn),-1,5.-eps) 

  rno = rnj 

  rno = where(rnj.lt.0.,rno,rnj) 
  rno = where(rnk.lt.0.,rno,rnk)  
  rno = where(rnl.lt.0.,rno,rnl)  
  rno = where(rnm.lt.0.,rno,rnm)  
  rno = where(rnn.lt.0.,rno,rnn)  

  slm = fild->slm

;;----------------------------------------------------------------------
;; new data structure for zonal average 
;;----------------------------------------------------------------------
 
 
  aa = rna(0,:,:) 
  bb = rnf(0,:,:) 
  cc = rni(0,:,:)  
 
  dd = rno(0,:,:)  

  prna = aa
  prnf = bb
  prni = cc

  prnj = dd
 
  pslm = slm(0,:,:) 

  ;prnb = where(pslm.gt.1.E-5,-999.,aa)    
  ;prng = where(pslm.gt.1.E-5,-999.,bb)    
  ;prnj = where(pslm.gt.1.E-5,-999.,cc)    
 

;;----------------------------------------------------------------------
;; open wks  
;;----------------------------------------------------------------------

  wks = gsn_open_wks(plotform,plotname)               ; open a ps file

;;----------------------------------------------------------------------
;; colormap
;;----------------------------------------------------------------------

  colormap = "testcmap"
  colormap = "amwg"
  colormap = "ViBlGrWhYeOrRe"
  colormap = "WhViBlGrYeOrRe"
  colormap2= "StepSeq25"
 
 ;gsn_define_colormap(wks,colormap)                   ; choose a colormap
  gsn_merge_colormaps(wks,colormap,colormap2)         ; choose a colormap
 ;gsn_draw_colormap(wks)
 ;exit

;;----------------------------------------------------------------------
;; plot settings:  
;;----------------------------------------------------------------------

  ;; for all panels 
 
  res                  = True    ; contour mods desired
  res@gsnFrame         = False   ; don't draw yet
  res@gsnDraw          = False   ; don't advance frame yet
  res@gsnMaximize      = False ;;True

  res@tiMainOn         = False   ; no title
  res@cnFillOn         = True    ; turn on color
  res@cnLinesOn        = False   ; turn off contour lines
  res@cnLineThicknessF = 1.     ; thicker lines
  res@cnInfoLabelOn    = False    ; no contour labels
  res@cnFillMode       = "RasterFill"  ; turn on raster mode
 ;res@cnFillDrawOrder  = "PreDraw"     ; draw contours before continents

  res@gsnMajorLonSpacing = 60.
  res@gsnMajorLatSpacing = 30.

  res@gsnSpreadColors = True    ; use full colormap

  if (colormap.eq."testcmap") then 
     res@gsnSpreadColorStart = 50           ; start at color 2
     res@gsnSpreadColorEnd   = 180          ; don't use added gray
  end if 
  if (colormap.eq."ViBlGrWhYeOrRe") then
     res@gsnSpreadColorStart = 24           ; start at color 2
     res@gsnSpreadColorEnd   = 95          ; don't use added gray
  end if
  if (colormap.eq."WhViBlGrYeOrRe") then
     res@gsnSpreadColorStart = 2            ; start at color 2
     res@gsnSpreadColorEnd   = 102          ; don't use added gray
  end if
 
  res@tiYAxisString  = ""  

;;----------------------------------------------------------------------
;; strings 
;;----------------------------------------------------------------------
  res@gsnLeftString    = "" 
  res@gsnCenterString  = "" 
  res@gsnRightString   = ""

  res@gsnLeftStringFontHeightF   = StringFontHeightF 
  res@gsnCenterStringFontHeightF = StringFontHeightF 
  res@gsnRightStringFontHeightF  = StringFontHeightF 

  res@gsnStringFont = gsnStringFont 
 
;;----------------------------------------------------------------------
;; lable bar 
;;----------------------------------------------------------------------
 
  res@lbLabelBarOn       = True       ; no individual label bar
  res@lbLabelStride      = 2          ; every other label bar label
  res@lbOrientation      = "Vertical" ; vertical label bar
  res@lbLabelFontHeightF = lbLabelFontHeightF 
  res@pmLabelBarWidthF   = 0.08        ; default is shorter
  res@pmLabelBarHeightF  = 0.5        ; default is taller

  res@trYReverse         = False      ; reverse y axis
 ;res@trYMinF            = 0.         ; set minimum Y-axis value
 ;res@trYMaxF            = 8000.      ; set maximum Y-axis value
 ;res@trXMinF            = 1949.      ; set minimum X-axis value
 ;res@trXMaxF            = 2006.      ; set maximum X-axis value

  res@vpWidthF           = 1.0        ; width of contour plots

;;----------------------------------------------------------------------
;; tickmark 
;;----------------------------------------------------------------------
;;res@tmYLMode = "Explicit"	
;;res@tmYLValues = (/100.,200.,400.,600.,800.,1000.,2000.,3000.,4000./)
;;res@tmYLLabels = (/"100","200","400","600","800","1000","2000","3000","4000"/)
;;res@tmYLMinorValues  = ispan(0,4000,200)

  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; axis  
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF

;;----------------------------------------------------------------------
;; plot array 
;;----------------------------------------------------------------------

  plot = new(np,graphic)                          ; create graphical array

;;----------------------------------------------------------------------
;; common setting 
;;----------------------------------------------------------------------

  res@cnLinesOn       = False    ; turn off contour lines
  res@cnFillOn        = True     ; turn on color
  res@cnLevelSelectionMode  = "ManualLevels"  ; set manual cn levels 
  res@cnMinLevelValF  = 0.2     ; set min contour level
  res@cnMaxLevelValF  = 2.0     ; set max contour level
  res@cnLevelSpacingF = 0.2
  res@lbLabelStride   =   1      ; every other label bar label  


;;----------------------------------------------------------------------
;; land panel 1 
;;----------------------------------------------------------------------

  ip = 0 

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prna,res) 
 
;;----------------------------------------------------------------------
;; land panel 2 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prnf,res) 
 
;;----------------------------------------------------------------------
;; land panel 3 
;;----------------------------------------------------------------------

  ip = ip + 1  

  res@gsnRightString  = gsnRightString(ip) 
  res@gsnLeftString   = gsnLeftString(ip) 
  res@gsnCenterString = gsnCenterString(ip) 

  plot(ip) = gsn_csm_contour_map_ce(wks,prni,res) 

;;----------------------------------------------------------------------
;; land panel 4
;;----------------------------------------------------------------------

  res@lbLabelFont      = "Helvetica-Bold"     ; label font
  res@lbLabelAlignment = "BoxCenters"         ; label orientation
  res@lbLabelPosition  = "Right"             ; label position
  res@lbLabelStrings =  (/"S3","Z","G","S2","S1"/) 

  ip = ip + 1

  res@gsnRightString  = "" ;;;gsnRightString(ip)
  res@gsnLeftString   = gsnLeftString(ip)
  res@gsnCenterString = gsnCenterString(ip)

  res@cnLevelSelectionMode  = "ExplicitLevels"  ; set manual cn levels 
  res@cnLevels              = (/ 1., 2., 3., 4. /) 
  res@gsnSpreadColors       = False 
  res@cnFillColors          = (/16, 107, 74, 5, 119/) 
  res@cnFillColors          = (/16, 107, 74, 6, 119/) 
  res@cnMonoFillPattern     = False
  res@cnFillPatterns        = (/9, 1, 10, 17, 8/)
  res@cnMonoFillScale      = False            ; want different densities 
  res@cnFillScales         = (/.1,.2,.3,.4,1,1,1.5,2.0,2.5/) ; the densities

  plot(ip) = gsn_csm_contour_map_ce(wks,prnj,res)
  
;;----------------------------------------------------------------------
;; panel plot 
;;----------------------------------------------------------------------

  resp                    = True   ; panel mods desired
  resp@gsnFrame           = False  ; don't advance frame yet
 ;resp@gsnPanelBottom     = 0.0    ; space for label bar
 ;resp@gsnPanelTop        = 0.5    ; only panel on lower half of page
 ;resp@gsnPanelLabelBar   = True   ; common label bar
 ;resp@pmLabelBarWidthF   = 0.8    ; label bar width
 
  resp@gsnPanelYWhiteSpacePercent = 5
 ;resp@gsnPanelXWhiteSpacePercent = 5
 
  gsn_panel(wks,plot(0:np-1),(/4,1/),resp)

  frame(wks) 

  system("ps2pdf "+plotfile) 
  ;system("mpack -s "+plotfile+" -a "+plotpdf+" k.zhang.iap@gmail.com") 

  ;; output plot information 


end
 
          


